#include <iostream>
#include <cmath>

#include "inv.h"

#define swap(a,b) {double temp=a; a=b; b=temp;}

std::vector<std::vector<double> > inverse(std::vector<std::vector<double> > a, unsigned n) {
    std::vector<std::vector<double> > inv(n, std::vector<double> (n));
    inv = a;
    std::vector<unsigned> indxc(n),indxr(n),ipiv(n);
    unsigned i,j,k,l,m,icol,irow;
    double big,dum,pivinv;

    irow=0;
    icol=0;

    for (j=0; j<n; j++) ipiv[j]=0;
    for (i=0; i<n; i++) {
        big=0;
        for (j=0; j<n; j++) {
            if (ipiv[j] != 1) {
                for (k=0; k<n; k++) {
                    if (ipiv[k] == 0) {
                        if (fabs(inv[j][k]) >= big) {
                            big = fabs(inv[j][k]);
                            irow=j;
                            icol=k;
                        }
                    } else if (ipiv[k] > 1) {
                        std::cout << "Сингулярная матрица 1!" << std::endl;
                        //a.clear();
                        //return a;
                    }
                }
            }
        }
        ++(ipiv[icol]);
        if (irow != icol) {
            for (l=0; l<n; l++) swap(inv[irow][l],inv[icol][l]);
        }

        indxr[i]=irow;
        indxc[i]=icol;

        if (inv[icol][icol] == 0) {
            std::cout << "Сингулярная матрица 2!" << std::endl;
            //a.clear();
            //return a;
        }

        pivinv=1/inv[icol][icol];
        inv[icol][icol]=1;

        for (l=0; l<n; l++) inv[icol][l] *= pivinv;
        for (m=0; m<n; m++) {
            if (m != icol) {
                dum = inv[m][icol];
                inv[m][icol]=0;
                for (l=0; l<n; l++) inv[m][l] -= inv[icol][l]*dum;
            }
        }
    }

    for (int p=n-1; p>=0; p--) {
        if (indxr[p] != indxc[p]) {
            for (k=0; k<n; k++) {
                swap(inv[k][indxr[p]],inv[k][indxc[p]])
            }
        }
    }

/*
    for (unsigned i=0; i<n; i++) {
        for (unsigned j=0; j<n; j++) {
            std::cout << "inv" << "[" << i << "]" << "[" << j << "] = " << inv[i][j] << std::endl;
        }
    }
*/
    return inv;
}
